o 

(N 

0^ 



Stochastic slowdown in evolutionary processes 

Philipp M. AltrockQ Chaitanya S. Gokhale, and Arne Traulsen 

Emmy-Noether Group for Evolutionary Dynamics, 

Department of Evolutionary Ecology, 

Max- Planck- Institute for Evolutionary Biology, 

August- Thienemann-Str. 2, D~24306 Plan, Germany 

(Dated: July 12, 2010) 

We examine birth-death processes with state dependent transition probabilities and at least one 
absorbing boundary. In evolution, this describes selection acting on two different types in a finite 
population where reproductive events occur successively. If the two types have equal fitness the 
system performs a random walk. If one type has a fitness advantage it is favored by selection, 
which introduces a bias (asymmetry) in the transition probabilities. How long does it take until 
advantageous mutants have invaded and taken over? Surprisingly, we find that the average time of 
such a process can increase, even if the mutant type always has a fitness advantage. We discuss 
this finding for the Moran process and develop a simplified model which allows a more intuitive 
understanding. We show that this effect can occur for weak but non-vanishing bias (selection) 
in the state dependent transition rates and infer the scaling with system size. We also address 
the Wright-Fisher model commonly used in population genetics, which shows that this stochastic 
slowdown is not restricted to birth-death processes. 

PACS numbers: 87.23.-n, 05.40.-a, 02.50.-r 
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INTRODUCION 



Birth-death processes belong to the most simple 
stochastic models and are applied in a variety of fields 
[iHB]. In physics these processes are connected e.g. to 
the study of one-dimensional classical diffusion in disor- 
dered media, anomalous transport, and molecular motors 
[TllTO]. In evolutionary biology, birth-death processes 
are commonly applied to model the evolution of traits 
with different reproductive fitness that are under natu- 
ral selection [21 |TT] . In the context of evolutionary game 
theory, this particular class of Markov chains has been 
used to model the spreading of successful strategies in a 
population of small size [T^20j . Naturally, the limit of 
weak selection is considered to be important in biology. 
It describes situations in which the effects of payoff dif- 
ferences are small, such that the evolutionary dynamics 
are mainly driven by random fluctuations. While this ap- 
proach has a long standing history in population genetics 
[2T| [22] , in the context of evolutionary game dynamics it 
has been introduced only recently tl4j. Often, from the 
discrete stochastic process a continuous limit or diffusion 
approximation is motivated, where typically the impact 
of the relevant parameters and time scales can be stud- 
ied more easily [TTll^M^ . Here, we consider the Moran 
process from theoretical population genetics and related 
processes. We address the speed of evolution when a res- 
ident population is taken over by mutants that are more 
fit. Under the low mutation rates that typically occur 
in biology, a mutant either goes extinct or takes over 
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the population before another mutation arises. Thus, for 
many purposes it is sufficient to address the evolution of 
two types in a one-dimensional system. 

In the following, we first recall general properties of 
birth-death processes (Sec. In]) and then address asym- 



metry in the transition probabilities (Sec. Ill I. In Sec. IV 



we then consider a more general Markov process to high- 
light that our main finding is not a special property of 
birth-death processes. 



II. STATE DEPENDENT BIRTH-DEATH 
PROCESS 



A one-dimensional birth-death process in position i 
can move to i — 1 or i -I- 1 with probabilities T[' and T^ . 
With probability 1 — T[' — T^ , the process stays in state 
i. We assume T^ = T^ — 0, such that i = and i = N 
are absorbing states. In discrete time, the probability to 
reach boundary N in t steps, starting from any i, obeys 
the master equation [6^ 
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The stationary conditional n"^ moment of Pf^ (t) is given 
by 
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The normalization constant, (/)f = J2't!Lo ^-f (^) ' ^^ ^^^ 
probability that the process gets absorbed at boundary 
N, called fixation probability in population genetics. For 
a recursion is obtained from Eq. M, (f)^ = (1 — 
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1^+0^1. With the boundary 
= 1, the solution reads [J] 



kN 



1 + V"^ n'' ^ 
1 + v^-^n'^ ^ 
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A measure for the duration of the process is the condi- 
tional mean time to absorption (average fixation time) 
r/^, i.e. the first moment of P.^{t). This gives the av- 
erage number of time steps until one of the two absorb- 
ing states is reached, starting from any i [3 [121 . A re- 
cursion for t/^ is obtained by multiplying each side of 
Eq. (fTl) with t and summing over all i [B], which yields 
c^frY^ (l-i;+~7^-)^frf + 7;-</.f_i(rfi + l)-f 
T'j (/'f^x('''i+i + !)• ^ similar recursion can be found for 
the conditional mean exit time t^, such that the mean 
life time of the process amounts to t° + rf . Solving 
recursively with the boundary conditions T(]^ = and 
T^ = 0, leads to the conditional mean time to reach 
state N, starting from i = 1, 
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One common example for a birth-death process with ab- 
sorbing states and N is the homogenous random walk, 
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< N and T^ 










This leads to 0f = i/N and rf = {N^ - l)/(6c). The 
reference case of population genetics is neutral evolution, 
where the symmetric transition probabilities are state de- 
pendent, T^ = i{N - i)/N'^. This resuhs in (pf = i/N 
andrf = N{N - 1) [51 [H]. 



the population. We consider the evolution of two types 
A and i? in a finite population of size N. Type A (with 
fitness Ja) is usually referred to as the mutant type, B 
(with fitness fs) is called the wild-type. Let i be the 
number of individuals of type A, such that iV — i is the 
number of B individuals. In general, the transition prob- 
abilities are 
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(5) 



In the following, we discuss different choices of /a and 
/b, as well as closely related, but simplified asymmetric 
transition rates. 



A. Constant fitness 

In the simplest case, the fitness of mutants is constant 
and does not depend on their abundance |llj . In our 
model, this can be parametrized as /^ = 1 -I- /? and fs — 
1 — p. In this case, the fixation probability of a single 
mutant is [TT] 



= (l-7)/(l-7") 



(6) 



where 7 = (1 — /3)/(l + /3). Up to linear order in /3 we 
have (/)f « N-'^ + I3{N - l)N-\ The larger the fitness 
advantage, the more likely the evolutionary takeover. For 
stronger selection (/3 > 0) an advantageous mutant is 
expected to fixate faster compared to neutral {[3 — 0). 



III. BIASED TRANSITION PROBABILITIES 

In this section, we examine how the state dependent 
transition probabilities influence the conditional mean 
exit time. We consider processes in which a parame- 
ter [3 continuously introduces a bias towards moving into 
one direction: For /3 — the transition probabilities are 
symmetric, T^ = T!^, but for /3 > 0, an asymmetry 
arises, T^ > T/^ . In evolutionary dynamics, /3 is usu- 
ally referred to as the intensity of selection. It governs 
the selective advantage (or disadvantage) of mutants in 
a wild-type population of finite size. Intuitively, it is 
clear that the time Ti does not depend trivially on /?, cf. 
Eq. Q. With increasing /?, the probability (j)f increases, 
but both 1/Tj^ and T~ jT'^ decrease in our setup. Thus, 
the average time t^ can increase or decrease with /?. In 
other words, despite increasing the tendency to move in 
the direction of a given boundary in each state, the con- 
ditional average time until this boundary is reached can 
still increase. 

In the Moran process, an individual selected for re- 
production proportional to fitness produces identical off- 
spring that replaces a randomly selected individual from 



B. Linear density dependence 

In general, the fitness of the two types will depend 
on their abundance. For example, the fitness / of each 
type can change linearly with i, Ja = I + l3{ai + b) 
and Jb = 1 — (3{ai + b). The bias (3 is bound such 
that fitness never becomes negative. Then, the transition 
probabilities are 
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N - l3{ai + h){N -2i) 
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(7) 



We have T^ 



^ N 



0, such that both boundaries are 
absorbing [HI [17]. For a < and aN + b>Q, type A is 
always fitter than type B, Ja > /b, but the conditional 
mean exit time t(^ is larger than neutral in a certain 
parameter range, compare FIG. IT] (a). In this case, a 
mutant that is fitter than the rest of the population needs 
more time to take over the population than a less fit 
mutant - intuitively, this should not be the case. The 
linear approximation of Tj'^ for /3 ^ N^^ (weak selection) 
reads 
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FIG. 1: (Color online) The conditional mean exit time 
Ti /ti ifi) (normalized) as a function of the bias (selection 
intensity) /?, or the mutation rate fi, for the four different 
models discussed in the main text. Symbols are simulations, 
lines show Eq. J4|. (a) Moran process with a = —0.1 and 
& = 2, see Eq. uV. (b) Parabolic-step process with i* = 11, 
Eq. (9|. (c) Constant -step process with i* = 9 and c — 0.5, 
Eq. I I2I). (d) Birt h-death process with directed mutations, 
Eqs. ( |15| and 1 16 1. The quantities f, P, and j3* indicate the 
maximal realtive increase of t^ , the according bias parame- 
ter, and the non-trivial value of /3 where t^ — rf (0), respec- 
tively (also compare FIG. [2|. The system size is TV = 20 in 
all panels, averages taken over 10^ realizations. 



see |5H1 m] . Note that the linear approximation of the 
conditional mean exit time depends only on the param- 
eter a, but not on 6, which holds for any system size. 
Hence, for small bias /3 and a < 0, the conditional aver- 
age time grows with increasing /3. This is an effect from 
state dependent fitness in finite populations, as it cannot 
occur for a = 0. 

The ratio T^ jT^ is a measure of the stochastic flow. 
Stochastic slowdown can occur if this ratio changes with 
the position (abundance of A) i, leading to an asymme- 
try. When /? becomes larger, t^ decreases again with /?, 
which is the strong selection behavior one would expect, 
compare FIG. fl] (a). 



C. Step— like asymmetry 

Is there a simpler process with similar characteristics? 
Indeed, we can introduce asymmetry also as a step in the 
fitness of the two types in our Moran process. This leads 
to parabolic transition probabilities with an additional 
step-like discontinuity, 






{\±^Q\i*-%\), 
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where Q\x\ is the step function (0[a; < 0] = and 0[a; > 
0] = 1). The integer i* is the location of the step. This 



process has the fixation probabilities 
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\ii> i* 



where i^\ 



(1 — 7)/(l — 7*^) is the probability to get from 
1 to fc, and 7 = (1 — /3)/(l + /3). Note that this general 
formula reduces to the standard fixation probability for 
constant fitness in the case of i* — N, cf. Eq.[6] For weak 
bias, /3 ^ l/N, we have 7 « 1 — 2/3, as well as 
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iii< i*, 
Hi > i*. 

(11) 



(f>f increases with f3 in this approximation, whereas 7 
decreases with /3. Hence, the mean exit time can also 
increase in an appropriate parameter range. The aver- 
age delay of the absorption is rather high in this case, 
cf. FIG. [1] (b), where it is 10%. FIG. [| (c) illustrates 
that even a delay of 400% is possible, but this delay de- 
creases with increasing i* . 

An even simpler model with stochastic slowdown is the 
constant -step process 



T^ 



c{l±l3e[i* ~i]) if 0< j<iV, 



(12) 



and T^ 






0, with i* < N, and the constant c cho- 
sen such that T^ + T~ < 1. Clearly, the fixation proba- 
bility of this process obeys Eqs. ( 10 ) and (111. Then, the 



remaining sums can be expressed by means of the exact 
form of (j)f , respecting that 1/T; only gives contribu- 
tions different from 1/c if ^ < i* . The conditional mean 
exit time r,^ can now be written in the form 
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With 7 ss 1 - 2/3 and Eq. ([llj) this leads to 

N^-1 {N -i*){N~l-i*)i*{l + i*) 
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/3. (14) 



The constant contribution is that of the homogenous ran- 
dom walk. The correction linear in /3 is always greater 
than or equal to zero, i.e. within the range of this approx- 
imation it just adds a positive value to the symmetric 
part. Also note that t{^ {P = 0,i*,c) serves as an upper 



bound for the mean exit time if i* > N — 1. Hence, be- 
low a certain threshold of the bias, t^^ is always greater 
than or equal to the homogenous random walk between 
absorbing boun darie s. This is surprising as the process 
defined by Eq. (fl2| fulfills T+ > T~ , and thus never 



gives a disadvantage to movement towards the bound- 
ary i = N. Moving into the direction of N is always 
at least as likely as moving into the opposite direction 
in this setup. In this particular process, the stochastic 
slowdown can be quite large, cf. FIG. [1] (c) and FIG. [2] 
(c). 

What is the effect of system size on this stochastic 
slowdown? Let /?* denote the upper bound of the pa- 
rameter f3 for which t^ [fi) > Ti {Q), which is the pa- 
rameter range in which slowdown can be observed. Ad- 
ditionally, with j3 we denote the parameter value of max- 
imal slowdown of the exit time t^ . They change with 
N and i* in both models with a step-like asymmetry, 
Eqs. ([9]) and ( |l2| ). The expansions linear in /3 are valid 
if 7V;9 < 1 [TsniTl W\- In FIG. [2] (a) and (b) we show 
that with increasing system size N , the quantities N/3{i*) 
and N/3*(i*) approach limiting curves if /3 is rescaled 
appropriately. Thus, stochastic slowdown does not rely 
on small system size, but /3* and /3 asymptotically scale 
as A'^"^. However, the maximal relative increase of the 
mean exit time itself, f — Ti {(i) / Ti {Q) , docs not scale 
with system size, f ^ iV", as illustrated in FIG.l2l(c). 



D. Directed mutations 

To stress the generality of the effect of stochastic slow- 
down in asymmetric birth-death processes we briefly dis- 
cuss a model with directed mutations. Fitness does not 
need to be position/state dependent to observe stochastic 
slowdown in population genetics. As above we consider 
two types, A and B, in a population of size A^, both hav- 
ing the same reproductive fitness. In one reproduction 
step of this Moran process, type B mutates to type A 
with a probability ^, back-mutations are excluded. This 
introduces asymmetry in the transition rates, 
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FIG. 2: (Color online) Scaling with system size for the 
two models with step like asymmetry: Parabolic-step model 
Eq. ||9| (FIG.[l](b)) on the left, constant-step model Eq. ([l2l 
with c = 1/2 (FIG. [1] (c)) on the right, (a) The threshold 
value iV/3*, defined by rf (/3*) = rf (0). Note that /3 < 1 per- 
mits a minimal value oii* /N only relatively far from zero, (b) 
A/3, defined as the bias parameter where the mean exit time 
tI' is maximal. When plotted against the asymmetry param- 
eter i* , both models approach a limit curve with growing size 
A. This suggests that non-trivial values of /3* and P can be 
found for any system size N after appropriate rescaling: The 
asymptotic scaling relations are /3 ~ N~^ , and /3* ~ A*'"^. 
(c) The maximal increase of the mean exit time (normal- 
ized), f = tI' {I3)/Tf (0), quickly approaches a limiting curve 
with growing N. This suggests the asymptotic scaling rela- 
tion f ~ A°. Open symbols A = 20, filled symbols A = 200, 
lines A = 2000. 



where i is the abundance of A. Obviously, T^ = Tq =0, 
but with directed mutations we have Tq*" > 0. The pro- 
cess has one absorbing boundary. The ratio of the tran- 
sition probabilities is T~/T^ w 1 — fiN/m, for mutation 
rates fi <C l/N'^. For larger /i, the dependence on the in- 
verse mutation rate makes the calculation of an approx- 
imation of Eq. (41) unwieldy. As fx increases we expect 
that A has an advantage during reproduction and hence, 
the conditional fixation time (that a single mutant takes 
over before going temporarily extinct) should decrease. 
Nevertheless, we observe an increase in the value of t^ , 



see FIG. R] (d). The time shows a maximum when /i is 
close to A^~^. 



A more general process is given in Appendix!^ There, 
we derive an expression for the fixation probability in a 
Wright -Fisher model with directed mutations. Although 
this quantity increases with /i, the associated conditional 
mean exit time also increases in a certain parameter 
range, compare FIG.|3] 



IV. STATE DEPENDENT WRIGHT-FISHER 
PROCESS 

The phenomenon of stochastic slowdown is not re- 
stricted to birth-death processes. It also occurs in the 
Wright-Fisher process that is commonly used in popula- 
tion genetics [TTJ [30] . Again, we consider a population of 
two types A and B. If i is the abundance of A, the fitness 
of each type is Ja — l + l3{ai + b), and /s — 1 — /3(ai + 6), 
respectively. Birth-death processes, such as the Moran 
model considered above, deal with one reproductive event 
at a time. Now, one time step of the Wright -Fisher pro- 
cess corresponds to one generation where all individuals 
reproduce: In each generation, the N individuals repro- 
duce a large number of offspring proportional to fitness. 
The new generation of size A^ is a random sample from 
this offspring pool, which corresponds to binomial sam- 
pling proportional to fitness. The transition probability 
to go from i to j A individuals reads |30j 



T'i^j — 



if A 



lfA + {N-l)fB 

jN - i) Ib 

tfA + {N- t)fE 
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(17) 



For this process, a closed treatment is not possible. Apart 
from simulations, for large N a diffusion approximation 
leads to analytical resuhs dH ISTHM] . With x = i/N, 
the process is approximately described by the Langevin 
equation dx = Di{x)dt + yjD2{x)dW{t), where W{t) is 
the Wiener process with zero mean and autocorrelation 
{W{t)W{s)) = min(t,s), [1]. The drift term Di{x) can 
be written as 



Di{x) = x{l-x)N- 



IaJx) - fsix) 
xfA{x) + {l-x)fB{x)' 

For the diffusion term D2{x) we find 

fA{x)fB{x) , Dl{x) 
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If the initial fraction of A types is Xg, the probability of 
absorption in a; = 1 (fixation probability) reads 



where 
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If there is no bias, /? = 0, we have /^(a;) = fsix) and 
hence Di{x) = 0. Thus, consistently with the previous 
section, we obtain (p{i/N) — i/N. For sufficiently weak 
bias, N/3 ^ 1, we have 



2Di{z) 
D2{z) 
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which leads to 

,. . , 2xo{l-Xo)N[aNil + xo)+3b] ^ 
(j){xo)^xo-\ p. (23) 

The conditional mean time this process takes to exit at 
a; = 1, t{xo), can be obtained from the associated back- 
ward Fokker-Planck equation [IT], 

Xo 1 
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For weak bias Eq. ( [22| ) holds, as well as S{x) 
2/3Nx^{aNx + 3b)l3. This results in 
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with 



Ci = a{7N^ + 13N + 6) + 186, 
C2 =6N{aN{N + 2) + 3b). 



For large A^, the right hand side of Eq. (26) simplifies, 
leading to 



r(l/A^) ^2N -l-a 



2N^{N ~3) 
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Hence, we can predict an increase of t(1/A^), in the case 
of state dependent bias with a < 0, also for the Wright- 
Fisher process, in particular when A always has a fit- 
ness advantage over B, see FIG.[3J This goes along with 
the findings for the Moran model in the previous sec- 
tion. Thus, the slowdown effect can also be observed in 
the traditional framework of population genetics, where 
times of fixation (or rather extinction) have been consid- 
ered typically for constant selection [TTl [55] . 



V. DISCUSSION 

This manuscript addresses several stochastic evolution- 
ary processes asking how long an advantageous muta- 
tion needs to take over. We have first concentrated on 
birth-death processes which model population dynam- 
ics with successive reproductive events, like the Moran 
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FIG. 3: (Color online) The conditional mean exit time (nor- 
malized) for the Wright-Fisher model with A^ — 1000, as 
a function of the rescaled bias (selection intensity, mutation 
rate). The line shows the analytical diffusion approximation 
result Eq. (24l, namely r(Ai'~^)/(2Ai'— 1). Symbols are simula- 
tion results. Left: The state dependent fitness model, Eq. ( 171 
(2x lO'' realizations, a — —0.1, b — N \a\). For relatively small 
bias /3 slow down is observed. Right: The directed mutations 
model, Eq. (All (5 x 10* realizations). Here, a strong slow- 



down effec can be observed over a wide range of the bias, 
Nn < 1. This is due to the different nature of the directed 
mutation process, which has only one absorbing boundary. 



Especially in the state dependent Moran or Wright- 
Fisher process, this can have a crucial impact on macro- 
scopic observable quantities. 
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process. However, the phenomenon of stochastic slow- 
down is also present in more general Markov processes, 
e.g. the Wright-Fisher process from population genet- 
ics. Stochastic slowdown is relevant in the invasion and 
fixation of beneficial traits with small state dependent 
selective advantage, which is typically assumed in evo- 
lutionary biology lijfjj. However, consequences of weak, 
but non- vanishing selection are hard to reveal in empir- 
ical studies, as the dynamics are still dominated by ran- 
dom genetic drift and averages over large ensembles are 
necessary. Biological examples of weak selection include 
amino acid substitutions which are only slightly advan- 
tageous or deleterious [57H55] . Weak state dependent 
fitness changes (such as the thresholds we discuss in our 
model with step-like asymmetry) may help explain situ- 
ations in which a substitution is likely, but takes a very 
long time. 

Our finding also has applications in evolutionary game 
theory ^D 42 : When a group of cooperative individuals 
is eventually driven to exinction by defectors, this process 
may take longer than the corresponding neutral process, 
although the defectors always have a fitness advantage. 
This observation is closely related to the fact that the 
conditional fixation time of an advantageous mutation is 
the same as the conditional fixation time of a deleterious 
mutation [551 133 . 

To sum up, we have shown that an asymmetric bias 
in a random walk, which is generic in population genet- 
ics, can lead to a counterintuitive observation that an 
advantageous mutant needs longer to take over the pop- 
ulation than a neutral mutant in the same system. This 
is a property of weakly biased systems, i.e. weak selec- 
tion, and is recovered for any system size if the intensity 
of selection is rescaled with A^"-*^. The relative maximal 
increase in time itself is independent of the system size. 



Appendix A: State dependent Wright— Fisher 
process with directed mutations 

Consider a finite population of size N, which consists 
of two types A and B. Both types have the same re- 
productive rate, which is set to one. In one generation, 
each type produces a large number of identical offspring 
proportional to its abundance. Additionally, a directed 
mutation from B to A can occur with probability /i. The 
next generation of size iV is a random sample from the 
offspring pool. The transition matrix reads 



T,- 



i 

N 



N 



N 



N - 



N 



-(1-m) 



N-j 



(Al) 



The conditional moments of this Markov chain are given 
by [IT] 



N 






(A2) 



In a diffusion approximation we rescale the state space 
as a; = i/N, and the timescale as At — 1/N, such that 
for large N the process is well described by the first two 
moments. 



N 



(A3) 



fc = 1, 2. For the given Markov chain Eq. ( Al ), the drift 
and diffusion terms read 



we obtain 



iiiy) ^ dz 



D2{z) 



lnD2{y) (A7) 



with 



Now, with 



/i((iV-l)/i + 2)-l- 



(A8) 



h{y) = cxp{-{h{y) - h{0))} = 



D2{v) 
D2iO)J ' 



we can calculate the second integral in Eq. (21 1 



S{x)^ / dyhiy) 



(A9) 



D^,+\x) - D!^+\0) 



D^(0)/i(2-^ + iV(2 + /i))-l' 



Hence, the fixation probability, Eq. (20), reads 



(f>{xo) 






(AlO) 



Di{x)^^iN{l~x), (A4) 

D2{x) = {l-x) ((1 -x)(N- l)n'^ + (1 - 2x)fi + x) . 

(A5) 

Next, we derive a closed expression for the probability 
that the process exits at a; = 1 without hitting the non- 
absorbing boundary x = first, starting from xq, 4>{xo), 
Eq. (20 1. The general expressions Eqs. (20), ([2T|), as well 



as Eqs. (^24L (25) hold. However, due to the different na- 



ture of this process, where only one absorbing boundary 
at a: = 1 exists, these quantities have a slightly different 
meaning. 

With 2Di{x)/D2{x) = 2Nfi/D2{x), where 



D2ix) = {l-x){N - l)n'^ + (1 - 2x)n + : 



(A6) 



As D2{0J = {{N - l)/i + 1)^, D2(l) = 1 - Ai, and 
Muifj^^Q D2{xo) = a;o, we have lim^^o0(a;o) = xq. Up to 
first order in mutation rate, we see that (/)(xo) increases 
with increasing bias, 

(f>{xo) Ki xo — (2iVa;o Inxo) /i. (All) 



With these expressions the conditional mean exit time, 
Eq. (24), can be tackled as well. However, we do not 



address the conditional mean exit time analytically, as 
its explicit form is elaborate and does not lead to fur- 
ther insight. From a numerical solution (Eq. ([24])) and 
from simulations the mean exit time of a single mutant, 
t{1/N), as a function of /j, is shown in FIG. [s] 
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